{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<h1 align=\"center\">Introduction to SimpleITKv4 Registration</h1>\n",
    "\n",
    "\n",
    "<table width=\"100%\">\n",
    "<tr style=\"background-color: red;\"><td><font color=\"white\">SimpleITK conventions:</font></td></tr>\n",
    "<tr><td>\n",
    "<ul>\n",
    "<li>Dimensionality and pixel type of registered images is required to be the same (2D/2D or 3D/3D).</li>\n",
    "<li>Supported pixel types are sitkFloat32 and sitkFloat64 (use the SimpleITK <a href=\"http://www.itk.org/SimpleITKDoxygen/html/namespaceitk_1_1simple.html#af8c9d7cc96a299a05890e9c3db911885\">Cast()</a> function if your image's pixel type is something else).\n",
    "</ul>\n",
    "</td></tr>\n",
    "</table>\n",
    "\n",
    "\n",
    "## Registration Components \n",
    "\n",
    "<img src=\"ITKv4RegistrationComponentsDiagram.svg\" style=\"width:700px\"/><br><br>\n",
    "\n",
    "There are many options for creating an instance of the registration framework, all of which are configured in SimpleITK via methods of the <a href=\"http://www.itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ImageRegistrationMethod.html\">ImageRegistrationMethod</a> class. This class encapsulates many of the components available in ITK for constructing a registration instance.\n",
    "\n",
    "Currently, the available choices from the following groups of ITK components are:\n",
    "\n",
    "### Optimizers\n",
    "\n",
    "The SimpleITK registration framework supports several optimizer types via the SetOptimizerAsX() methods, these include:\n",
    "\n",
    "<ul>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ExhaustiveOptimizerv4.html\">Exhaustive</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1AmoebaOptimizerv4.html\">Nelder-Mead downhill simplex</a>, a.k.a. Amoeba.\n",
    "  </li>\n",
    "  <li>\n",
    "   <a href=\"https://itk.org/Doxygen/html/classitk_1_1PowellOptimizerv4.html\">Powell optimizer</a>.\n",
    "  </li>\n",
    "  <li>\n",
    "   <a href=\"https://itk.org/Doxygen/html/classitk_1_1OnePlusOneEvolutionaryOptimizerv4.html\">1+1 evolutionary optimizer</a>.\n",
    "  </li>\n",
    "  <li>\n",
    "  Variations on gradient descent:\n",
    "  <ul>\n",
    "    <li>\n",
    "    <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1GradientDescentOptimizerv4Template.html\">GradientDescent</a>\n",
    "    </li>\n",
    "    <li>\n",
    "    <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1GradientDescentLineSearchOptimizerv4Template.html\">GradientDescentLineSearch</a>\n",
    "    </li>\n",
    "    <li>\n",
    "    <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1RegularStepGradientDescentOptimizerv4.html\">RegularStepGradientDescent</a>\n",
    "    </li>\n",
    "  </ul>\n",
    "  </li>\n",
    "  <li>\n",
    "    <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ConjugateGradientLineSearchOptimizerv4Template.html\">ConjugateGradientLineSearch</a> \n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1LBFGSBOptimizerv4.html\">L-BFGS-B</a> (Limited memory Broyden,  Fletcher,Goldfarb,Shannon-Bound Constrained) - supports the use of simple constraints ($l\\leq x \\leq u$)  \n",
    "  </li>\n",
    "</ul>\n",
    "\n",
    " \n",
    "### Similarity metrics\n",
    "\n",
    "The SimpleITK registration framework supports several metric types via the SetMetricAsX() methods, these include:\n",
    "\n",
    "<ul>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1MeanSquaresImageToImageMetricv4.html\">MeanSquares</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1DemonsImageToImageMetricv4.html\">Demons</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1CorrelationImageToImageMetricv4.html\">Correlation</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1ANTSNeighborhoodCorrelationImageToImageMetricv4.html\">ANTSNeighborhoodCorrelation</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1JointHistogramMutualInformationImageToImageMetricv4.html\">JointHistogramMutualInformation</a>\n",
    "  </li>\n",
    "  <li>\n",
    "  <a href=\"http://www.itk.org/Doxygen/html/classitk_1_1MattesMutualInformationImageToImageMetricv4.html\">MattesMutualInformation</a>\n",
    "  </li>\n",
    "</ul>\n",
    "\n",
    "\n",
    "### Interpolators\n",
    "\n",
    "The SimpleITK registration framework supports several interpolators via the SetInterpolator() method, which receives one of\n",
    "the <a href=\"http://www.itk.org/SimpleITKDoxygen/html/namespaceitk_1_1simple.html#a7cb1ef8bd02c669c02ea2f9f5aa374e5\">following enumerations</a>:\n",
    "<ul>\n",
    "<li> sitkNearestNeighbor </li>\n",
    "<li> sitkLinear </li>\n",
    "<li> sitkBSpline </li>\n",
    "<li> sitkGaussian </li>\n",
    "<li> sitkHammingWindowedSinc </li>\n",
    "<li> sitkCosineWindowedSinc </li>\n",
    "<li> sitkWelchWindowedSinc </li>\n",
    "<li> sitkLanczosWindowedSinc </li>\n",
    "<li> sitkBlackmanWindowedSinc </li>\n",
    "</ul>\n",
    "\n",
    "## Data -  Retrospective Image Registration Evaluation\n",
    "\n",
    "We will be using part of the training data from the Retrospective Image Registration Evaluation (<a href=\"http://www.insight-journal.org/rire/\">RIRE</a>) project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "library(SimpleITK)\n",
    "library(ggplot2)\n",
    "# Utility method that either downloads data from the Girder repository or\n",
    "# if already downloaded returns the file name for reading from disk (cached data).\n",
    "source(\"downloaddata.R\")\n",
    "\n",
    "# Always write output to a separate directory, we don't want to pollute the source directory. \n",
    "OUTPUT_DIR = \"Output\"\n",
    "\n",
    "# alpha blending value we will use for displaying an overlay of the registered images.\n",
    "alpha <- 0.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility functions\n",
    "A number of utility callback functions for plotting the similarity metric during registration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Callback invoked when the StartEvent happens, sets up our new data.\n",
    "# Functions ending in _jn are for use with Jupyter notebooks, as the display\n",
    "# behaviour is a bit different.\n",
    "# Note that we could use a special plotting environment instead of the global\n",
    "# environment, but we may want to use the metrics elsewhere and they are easier to\n",
    "# access if they are global\n",
    "start_plot <- function()\n",
    "{   #global empty vectors (via assignment operator) \n",
    "    metric_values <<- c()\n",
    "    multires_iterations <<- c()\n",
    "}\n",
    "\n",
    "# Callback invoked when the EndEvent happens, do cleanup of global data.\n",
    "end_plot <- function()\n",
    "{    \n",
    "    rm(metric_values, pos = \".GlobalEnv\")\n",
    "    rm(multires_iterations, pos = \".GlobalEnv\")\n",
    "}\n",
    "# In notebooks we only get one plot per cell, so use the end to plot it\n",
    "end_plot_jn <- function()\n",
    "{    \n",
    "    Multi <- rep(NA, length(metric_values))\n",
    "    Multi[multires_iterations] <- \"M\"\n",
    "    DDF <- data.frame(IterationNumber=1:length(metric_values), \n",
    "                      MetricValue=metric_values,\n",
    "                      MultiresIteration=Multi)\n",
    "    DDFM <- subset(DDF, !is.na(MultiresIteration))\n",
    "    pl <- ggplot(DDF, aes(x=IterationNumber, y=MetricValue)) + \n",
    "          geom_line() + \n",
    "          geom_point(data=DDFM, aes(colour=MultiresIteration)) +\n",
    "          theme(legend.position=\"none\")    \n",
    "    print(pl)\n",
    "    rm(metric_values, pos = \".GlobalEnv\")\n",
    "    rm(multires_iterations, pos = \".GlobalEnv\")\n",
    "}\n",
    "# Callback invoked when the IterationEvent happens, update our data and display new figure.\n",
    "# Note that this won't appear as an animation in R studio, but you can use the arrows to cycle\n",
    "# through plots\n",
    "plot_values <- function(registration_method)\n",
    "{\n",
    "    metric_values <<- c(metric_values, registration_method$GetMetricValue())\n",
    "    Multi <- rep(NA, length(metric_values))\n",
    "    Multi[multires_iterations] <- \"M\"\n",
    "    DDF <- data.frame(IterationNumber=1:length(metric_values),\n",
    "                      MetricValue=metric_values,\n",
    "                      MultiresIteration=Multi)\n",
    "    DDFM <- subset(DDF, !is.na(MultiresIteration))\n",
    "\n",
    "    pl <- ggplot(DDF, aes(x=IterationNumber, y=MetricValue)) +\n",
    "          geom_line() +\n",
    "          theme(legend.position=\"none\")\n",
    "\n",
    "    if(nrow(DDFM) > 1) {\n",
    "        pl <- pl + geom_point(data=DDFM, aes(colour=MultiresIteration))\n",
    "    }\n",
    "    print(pl)\n",
    "    dev.flush()\n",
    "    Sys.sleep(0)\n",
    "}\n",
    "# Use this one inside a notebook\n",
    "plot_values_jn <- function(registration_method)\n",
    "{    \n",
    "    ## No point attempting to plot every one in a notebook\n",
    "    metric_values <<- c(metric_values, registration_method$GetMetricValue())\n",
    "}\n",
    "                                 \n",
    "# Callback invoked when the sitkMultiResolutionIterationEvent happens, update the index into the \n",
    "# metric_values list. \n",
    "update_multires_iterations <- function()\n",
    "{\n",
    "    multires_iterations <<- c(multires_iterations, length(metric_values)+1)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read images\n",
    "\n",
    "We first read the images, casting the pixel type to that required for registration (Float32 or Float64) and look at them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "fixed_image <- ReadImage(fetch_data(\"training_001_ct.mha\"), \"sitkFloat32\")\n",
    "moving_image <- ReadImage(fetch_data(\"training_001_mr_T1.mha\"), \"sitkFloat32\")\n",
    "\n",
    "Show(fixed_image,\"fixed_image\")\n",
    "Show(moving_image, \"moving_image\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Initial Alignment\n",
    "\n",
    "Use the CenteredTransformInitializer to align the centers of the two volumes and set the center of rotation to the center of the fixed image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "initial_transform <- CenteredTransformInitializer(fixed_image, \n",
    "                                                  moving_image, \n",
    "                                                  Euler3DTransform(), \n",
    "                                                  \"GEOMETRY\")\n",
    "\n",
    "moving_resampled <- Resample(moving_image, fixed_image, initial_transform)\n",
    "Show((1-alpha)*fixed_image + alpha*moving_resampled, \"initial alignment\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Registration\n",
    "\n",
    "The specific registration task at hand estimates a 3D rigid transformation between images of different modalities. There are multiple components from each group (optimizers, similarity metrics, interpolators) that are appropriate for the task. Note that each component selection requires setting some parameter values. We have made the following choices:\n",
    "\n",
    "<ul>\n",
    "<li>Similarity metric, mutual information (Mattes MI):\n",
    "<ul>\n",
    "  <li>Number of histogram bins, 50.</li>\n",
    "  <li>Sampling strategy, random.</li>\n",
    "  <li>Sampling percentage, 1%.</li>\n",
    "</ul>\n",
    "</li>\n",
    "<li>Interpolator, sitkLinear.</li>\n",
    "<li>Optimizer, gradient descent: \n",
    "<ul>\n",
    "  <li>Learning rate, step size along traversal direction in parameter space, 1.0 .</li>\n",
    "  <li>Number of iterations, maximal number of iterations, 100.</li>\n",
    "  <li>Convergence minimum value, value used for convergence checking in conjunction with the energy profile of the similarity metric that is estimated in the given window size, 1e-6.</li>\n",
    "  <li>Convergence window size, number of values of the similarity metric which are used to estimate the energy profile of the similarity metric, 10.</li>\n",
    "</ul>\n",
    "</li>\n",
    "</ul>\n",
    "\n",
    "\n",
    "Perform registration using the settings given above, and take advantage of the built in multi-resolution framework, use a three tier pyramid.  \n",
    "\n",
    "In this example we plot the similarity metric's value during registration. Note that the change of scales in the multi-resolution framework is readily visible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "registration_method <- ImageRegistrationMethod()\n",
    "\n",
    "# Similarity metric settings.\n",
    "registration_method$SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)\n",
    "registration_method$SetMetricSamplingStrategy(\"RANDOM\")\n",
    "registration_method$SetMetricSamplingPercentage(0.01)\n",
    "\n",
    "registration_method$SetInterpolator(\"sitkLinear\")\n",
    "\n",
    "# Optimizer settings.\n",
    "registration_method$SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, convergenceMinimumValue=1e-6, convergenceWindowSize=10)\n",
    "registration_method$SetOptimizerScalesFromPhysicalShift()\n",
    "\n",
    "# Setup for the multi-resolution framework.            \n",
    "registration_method$SetShrinkFactorsPerLevel(shrinkFactors = c(4,2,1))\n",
    "registration_method$SetSmoothingSigmasPerLevel(smoothingSigmas=c(2,1,0))\n",
    "registration_method$SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()\n",
    "\n",
    "# Don't optimize in-place, we would possibly like to run this cell multiple times.\n",
    "registration_method$SetInitialTransform(initial_transform, inPlace=FALSE)\n",
    "\n",
    "# Connect all of the observers so that we can perform plotting during registration. \n",
    "# (assignment of the return values to dev_null so they aren't printed)\n",
    "dev_null <- registration_method$AddCommand(\"sitkStartEvent\", start_plot)\n",
    "dev_null <- registration_method$AddCommand(\"sitkEndEvent\", end_plot_jn)\n",
    "dev_null <- registration_method$AddCommand(\"sitkMultiResolutionIterationEvent\", update_multires_iterations) \n",
    "dev_null <- registration_method$AddCommand(\"sitkIterationEvent\", function() plot_values_jn(registration_method))\n",
    "\n",
    "final_transform <- registration_method$Execute(fixed_image, moving_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post registration analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Query the registration method to see the metric value and the reason the optimization terminated. \n",
    "\n",
    "The metric value allows us to compare multiple registration runs as there is a probabilistic aspect to our registration, we are using random sampling to estimate the similarity metric.\n",
    "\n",
    "Always remember to query why the optimizer terminated. This will help you understand whether termination is too early, either due to thresholds being too tight, early termination due to small number of iterations - numberOfIterations, or too loose, early termination due to large value for minimal change in similarity measure - convergenceMinimumValue)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cat(paste0(\"Final metric value: \",registration_method$GetMetricValue(),\"\\n\"))\n",
    "cat(paste0(\"Optimizer\\'s stopping condition, \",registration_method$GetOptimizerStopConditionDescription()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now visually inspect the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "simpleitk_error_allowed": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "moving_resampled <- Resample(moving_image, fixed_image, final_transform)\n",
    "Show((1-alpha)*fixed_image + alpha*moving_resampled, \"final alignment\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are satisfied with the results, save them to file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "WriteImage(moving_resampled, file.path(OUTPUT_DIR, \"RIRE_training_001_mr_T1_resampled.mha\"))\n",
    "WriteTransform(final_transform, file.path(OUTPUT_DIR, \"RIRE_training_001_CT_2_mr_T1.tfm\"))"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.2.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
